NASA/TM-2005-2 1 3764 



Factorization of the Compressible Navier- 
Stokes Equations 


Thomas W. Roberts 

Langley Research Center, Hampton, Virginia 


December 2005 



The NASA STI Program Office ... in Profile 


Since its founding, NASA has been dedicated to the 
advancement of aeronautics and space science. The 
NASA Scientific and Technical Information (STI) 
Program Office plays a key part in helping NASA 
maintain this important role. 

The NASA STI Program Office is operated by 
Langley Research Center, the lead center for NASA’s 
scientific and technical information. The NASA STI 
Program Office provides access to the NASA STI 
Database, the largest collection of aeronautical and 
space science STI in the world. The Program Office is 
also NASA’s institutional mechanism for 
disseminating the results of its research and 
development activities. These results are published by 
NASA in the NASA STI Report Series, which 
includes the following report types: 

• TECHNICAL PUBLICATION. Reports of 
completed research or a major significant phase 
of research that present the results of NASA 
programs and include extensive data or 
theoretical analysis. Includes compilations of 
significant scientific and technical data and 
information deemed to be of continuing 
reference value. NASA counterpart of peer- 
reviewed formal professional papers, but having 
less stringent limitations on manuscript length 
and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 
and technical findings that are preliminary or of 
specialized interest, e.g., quick release reports, 
working papers, and bibliographies that contain 
minimal annotation. Does not contain extensive 
analysis. 

• CONTRACTOR REPORT. Scientific and 
technical findings by NASA-sponsored 
contractors and grantees. 


• CONFERENCE PUBLICATION. Collected 
papers from scientific and technical 
conferences, symposia, seminars, or other 
meetings sponsored or co-sponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 
technical, or historical information from NASA 
programs, projects, and missions, often 
concerned with subjects having substantial 
public interest. 

• TECHNICAL TRANSLATION. English- 
language translations of foreign scientific and 
technical material pertinent to NASA’s mission. 

Specialized services that complement the STI 
Program Office’s diverse offerings include creating 
custom thesauri, building customized databases, 
organizing and publishing research results ... even 
providing videos. 

For more information about the NASA STI Program 
Office, see the following: 

• Access the NASA STI Program Home Page at 
http://www. sti. nasa.gov 

• E-mail your question via the Internet to 
help@sti.nasa.gov 

• Fax your question to the NASA STI Help Desk 
at (301) 621-0134 

• Phone the NASA STI Help Desk at 
(301) 621-0390 

• Write to: 

NASA STI Help Desk 

NASA Center for AeroSpace Information 

7121 Standard Drive 

Hanover, MD 21076-1320 



NASA/TM-2005-2 1 3764 



Factorization of the Compressible Navier- 
Stokes Equations 


Thomas W. Roberts 

Langley Research Center, Hampton, Virginia 


National Aeronautics and 
Space Administration 

Langley Research Center 
Hampton, Virginia 23681-2199 


December 2005 



Available from: 


NASA Center for AeroSpace Information (CASI) 
7121 Standard Drive 
Hanover, MD 21076-1320 
(301) 621-0390 

National Technical Information Service (NTIS) 
5285 Port Royal Road 
Springfield, VA 22161-2171 
(703) 605-6000 



Contents 

Abstract iv 

1 Introduction 1 

2 The Navier-Stokes equations 2 

3 Relaxation and principal linearization 2 

3.1 Principal linearization of inviscid terms 3 

3.2 Principal linearization of the viscous terms 4 

3.3 The complete principal linearization 5 

3.4 Interpretation of the principal linearization 6 

4 Factors of the Navier-Stokes equations 6 

4.1 Inviscid equations 7 

4.2 Low speed flow 8 

4.3 A special case of the Prandtl number 9 

5 Modes of the principal linearization 9 

6 Conclusions 11 


iii 



Abstract 


The Navier-Stokes equations for a Newtonian ideal gas are examined to determine the factorizable form of 
the equations relevant to the construction of a factorizable relaxation scheme. The principal linearization 
of the equations is found by examining the relative magnitude of the terms for short-wavelength errors. 
The principal part of the operator is then found. Comparison of the factors of the Navier-Stokes and Euler 
equations differ qualitatively because of the coupling of entropy and pressure through thermal diffusion. 
Special cases of the factorization are considered. 


IV 



1 Introduction 


Considerable progress has been made in recent years in the development of fast multigrid methods for the 
Euler and Navier-Stokes equations. In particular, the development of factorizable discretizations shows the 
promise of attaining ideal, or textbook, multigrid efficiency. A factorizable discretization is one in which the 
governing equations are partitioned into component elliptic, hyperbolic and parabolic subsystems. Each of 
these subsystems is then solved or relaxed by an appropriate numerical algorithm. By treating each partition 
optimally, the full system can be solved with the efficiency of the slowest converging subsystem. These ideas 
have been propounded by Achi Brandt for some time (Brandt, 1985) and are the subject of a recent review 
by Thomas, Diskin, and Brandt (2003). 

In the multigrid literature, factorizability is described in terms of the properties of the determinant of 
the system of equations (Sidilkover, 1999; Ta’asan, 1993; Thomas et ah, 2003). To relax the system of 
equations, the increment in the dependent variables is written in terms of a set of auxiliary variables, or 
“ghost” variables in Brandt’s terminology. The transformation of variables is designed such that the system 
of equations in terms of the auxiliary variables is in triangular form, with the diagonal elements being the 
factors of the determinant. Increments in the auxiliary variables are used to update the original variables 
though the transformation. The triangularization of the system by means of the ghost variables is nothing 
more than decomposing the system into a set of nearly linearly independent modes. The factors of the 
determinant are the operators which act on each mode, and the auxiliary variables are the mode amplitudes. 
The key point is that the modes represent the short-wavelength components of the solution error. Only the 
terms in the equations that contribute strongly at the short wavelengths need to be considered in factoring 
the equations. These terms are called the principal part of the operator. 

The factorizable property may be examined in a variety of ways. Thomas et al. (2003) view the con- 
struction of a factorizable discrete operator by treating the relaxation process as being an approximation to 
a Newton iteration. Starting with a linearization of the full system, they eliminate the subprincipal terms 
of the Jacobian and express the solution updates in terms of ghost variables. Ta’asan (1993) works directly 
with the full nonlinear Euler equations to derive a set of “canonical variables” corresponding to the different 
partitions of the equations. Using these variables he constructs a multigrid solver based on direct relaxation 
of the nonlinear system (Ta’asan, 1994). Sidilkover (1999) devises a factorizable scheme for the Euler equa- 
tions by designing the difference operators such that the discrete potential, entropy and vorticity modes are 
preserved exactly in the constant coefficient case. This is generalized to a conservative discretization of the 
full nonlinear system in a straightforward way (Roberts, Sidilkover, and Thomas, 2000). As with Ta’asan, 
Sidilkover prefers to directly relax the full nonlinear system. 

Although factorizable schemes have been constructed and demonstrated textbook multigrid efficiency for 
the incompressible Navier-Stokes equations (Brandt & Yavneh, 1992; Thomas, Diskin, and Brandt, 2001; 
Swanson, 2001) and for the incompressible and compressible Euler equations (Roberts et al., 2000; Roberts, 
Swanson, and Sidilkover, 1999; Ta’asan, 1994) in two dimensions, the development of such schemes for the 
compressible Navier-Stokes equations is still in its infancy. An initial attempt at such a scheme has been 
reported by Thomas, Diskin, and Brandt (1999). A difficulty arises because the nature of the factors in 
the viscous compressible case is qualitatively different than either the viscous incompressible or the inviscid 
compressible equations. Thermal conductivity has the effect of coupling the energy and continuity equations 
in compressible viscous flow. For certain values of the Prandtl number, described below, the two factors 
can be decoupled, and in general the coupling is weak. Thomas et al. (2003) present one possible approach 
to factorizing the compressible Navier-Stokes equations. Relaxation of their system will require a block 
relaxation of the two thermodynamic variables, which are pressure and internal energy in their case. A 
possible disadvantage of their formulation is that the energy and continuity modes do not decouple even for 
the appropriate value of the Prandtl number. 

In the remainder of this paper, the principal linearization of the compressible Navier-Stokes is developed. 
The full Navier-Stokes equations are presented in Section 2. Following Thomas et al. (2003), in Section 3 
the motivation is given by viewing relaxation methods as approximate Newton solvers. The principal lin- 
earization is derived from a perturbation expansion of the full equations. A discussion of various limits and 
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special cases of the factorization is presented in Section 4. The inviscid and incompressible limits are shown 
to lead to the appropriate factorization of the Euler and incompressible Navier-Stokes equation. With this 
factorization in mind, a decomposition into vortical, entropy and acoustic modes is developed, which may 
form the basis of a fast multigrid solver, in Section 5. These modes transform smoothly to the inviscid case 
as well. In addition, this decomposition has a counterpart in the study of compressible turbulent flow. 


2 The Navier-Stokes equations 


The analysis which follows is for a Newtonian fluid obeying an ideal gas law, with constant specific heats. 
Only steady flows are considered. For the purposes of examining the factors of the Navier-Stokes equations, 
the primitive variables are the most convenient. Letting (■ Ui,p,s ) be the state vector of velocity, pressure, 
and entropy, the Navier-Stokes equations may be written 


1 

UjUij + -p^ 

P 

pc 2 Ujj + UjP.j 


1—1 | 

(la) 

(7 - 1) ($ - q jt j ) , 

(lb) 

^ - Qjj) > 

(lc) 


where p is the density, T is the static temperature and 7 is the ratio of the specific heats. The Einstein 
summation convention is used, and partial differentiation with respect to the i-th coordinate is indicated by 
the comma subscript notation. The viscous stress tensor is defined as Ty = p'dkkdij + 2 pdij, where Sij is the 
Kronecker delta and di 3 is the deviatoric strain tensor d l:] = (u-- h] + Ujj ) /2. The first and second coefficients 
of viscosity p and p' are related through the definition of the bulk viscosity, 3 k = 3 p' + 2p > 0. By Stokes’ 
hypothesis k = 0, which gives p' = —2p/2>. However, it is not necessary to use this relation to determine the 
principal linearization of the equations. 

The two terms on the right-hand side of the entropy equation contain the effects of viscous heating and 
thermal diffusion. The dissipation function <I> = d l:j T t j > 0 is the viscous work done on the fluid. The heat 
flux vector is given by qi = —AT* = — (pc p /Pv) T*, where A is the thermal conductivity, c p is the specific 
heat at constant pressure, and Pr is the Prandtl number. The temperature and density may be expressed 
in terms of the primary thermodynamic variables s and p: 


T 

T~o 

P_ 
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(7-1 )h 

g(s-so )/c p 


/ \ 1/7 

[ P_ ] e -(s-so)/cp 

\PoJ 


where the subscript 0 denotes a reference value. The equation of state can also be written p = ((7 — 1) /j) pc p T , 
and the speed of sound is given by c 2 = 7 p/ p = (7 — 1) c p T. 


3 Relaxation and principal linearization 

Relaxation methods are a class of iterative methods for the solution of a system of algebraic equations. An 
iterative method for solving a system of equations Aip = b consists of generating a sequence of approximations 
{tpo VA P 2 ■ ■ •} such that liuin^oo ip n = (p. The system is linearized about the n-th iterate to obtain an 
equation for the error, A v ((p — <p n ) mb — A(p n where A v = d (Aip)/dip is the Jacobian of A and b — Ap n is 
the residual at the n-th iteration. When the system of equations is linear, then A v = A 1 and the equation 
for the error is exact. 

A Newton solver is obtained when all the terms of the linearization are kept and the resulting system 
for the error is solved exactly. For discretizations of partial differential equations, such as the Navier-Stokes 
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equations, this requires the inversion of a sparse, banded matrix. Newton methods have the outstanding 
property of quadratic convergence once the approximate solution is near the exact solution. Unfortunately, 
direct inversion is expensive and often impractical. Furthermore, the initial iterates of Newton’s method can 
be far from the exact solution — for nonlinear problems, many iterations may be required before the region 
of quadratic convergence is reached. 

Relaxation methods avoid direct inversion by replacing the Jacobian A v with an approximation that 
is inexpensive to invert. Generally, this means throwing away all or some of the off-diagonal terms of the 
Jacobian. For example, keeping just the diagonal terms leads to Jacobi iteration; the Gauss-Seidel method 
results from throwing away the terms above the diagonal. Inverting the matrix involves only the cost of a 
matrix multiplication, which because of the sparsity of the matrix is 0(n). 

Because it involves only local operations, relaxation is effective in eliminating the short-wavelength com- 
ponents of the error, where the short wavelengths are those on the order of h, the grid spacing. For this reason 
only those terms of the linearization that contribute strongly to the short- wavelength error components need 
be kept. These terms are referred to as the “principal part” of the linearization. Efficient multigrid methods 
are based on a relaxation method, or smoother, for reducing the short-wavelength component of the error, 
and a coarse-grid operator to reduce the long- wavelength error. 

In Thomas et al. (2003) the authors describe the process of determining the principal part by reference 
to the operator alone. Taking the determinant, only those terms that contribute to the leading order are 
retained. In the following treatment, the principal part is found by expanding the equations in a small 
perturbation about a smooth state and keeping the dominant terms. This leads to the same result as 
Thomas, et ah, but has the advantage of showing directly why certain terms are subprincipal. Furthermore, 
the procedure followed here is more familiar to fluid dynamicists and is very similar to a perturbation theory 
used in the study of compressible turbulence (Chu & Kovasznay, 1958; Kovasznay, 1953). 

3.1 Principal linearization of inviscid terms 

The principal linearization of the Navier-Stokes equations starts by considering a full linearization and 
keeping only those terms that are dominant on the scale of the grid spacing. The linearization is found by 
replacing the variables (s,w.j,p) with 


Ui + e-Ui, 

(2a) 

p + ep, 

(2b) 

s + es 

(2c) 


where e< 1. Equations (1) are expanded in e and the linear terms are kept. 

The inviscid terms of the Navier-Stokes equations (1) are those terms on the left-hand side. Let U be the 
velocity scale and L be the physical length scale of the flow. These are local scales, and can be different in 
different regions of the flow. For example, L can be a chord length, leading-edge radius, or boundary-layer 
thickness. The grid spacing is h. It is assumed that /i«L, or in other words, the grid is sufficiently fine to 
resolve the physical features of the flow. Making the substitutions Eqs. (2) into Eqs. (1) and keeping only 
the term of O(e) yields the linearization of the inviscid terms 


, 1 ~ P P,i , s P,i 

UjUij + UjUij + -Pi 1 , 

(3a) 

p pc 2 P Cp p 


pC 2 Uj } j + 'ypUjj + UjPj +PjUj, 

(3b) 

UjS,j + UjSj. 

(3c) 


The scaling of the variables is 

Ui,Ui~U p,p ~ pU 2 s,s~c p M 2 . 
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Because only the short-wavelength components of the error are of interest, the derivatives of the perturbation 
quantities scale as 


UiJ 




c p M 2 

h 


but the derivatives of the solution values scale as 


U 

L 



c p M 2 
L 


Using the assumption that h -C L, it is seen that terms involving the derivatives of solution quantities 
can be neglected compared to the terms involving derivatives of the perturbation quantities. This gives the 
principal linearization of the inviscid terms 


1 „ 

+ - p,i , 
p 

(4a) 


(4b) 

Ai- 

(4c) 


Also, note that the ratio of the second term to the first term of Eq. (4b) is 0{M 2 ). Thus for low Mach 
number flows, M 2 <C 1, the pressure advection term in the continuity equation is subprincipal. 


3.2 Principal linearization of the viscous terms 

The right-hand sides of Eqs. (lc) and (lb), the entropy and continuity equations, consist of the heat diffusion 
and viscous dissipation terms. The principal part of the linearization of each term is considered. The scalings 
of the viscosity and temperature are 


M, p' 


pUL 

Re 


T,T ~ 


U2 

Cp 


The Reynolds number Re is based on the relevant local velocity and length scales. Also, it is assumed 
that dp/dT p/T, d 2 p/dT 2 <C p/T 2 , so the variations of the viscosity can be neglected when linearizing. 
Examining the scaling of the deviatoric strain tensor and its derivatives gives 


dij 







U 

h 2 ' 


As a consequence, the dominant term in the linearization of the dissipation function $ is 


2p'dudjj ■ \pd.jdij. 


Considering the thermal diffusion term next, note that the temperature derivatives scale as 
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from which it follows that the dominant term in the heat-flux linearization is 

Pr ' ' 

Using the equation of state to express this in terms of s and p gives 


v / pc “ Sy 


Pr \7 - 1 Cp 


■P,i 
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Examination of the relative magnitude of the two terms yields 


2 n' dudjj +4 fidijdij 


P£pf 

Pr ’ 


pU 3 L 1 
Re Lh ’ 
pU 3 L 1 
Re Pr E? 


It is seen that that thermal diffusion term is larger by 0{Pr~ 1 L/h), so the dissipation term is subprincipal 
and can be neglected in the linearization. 

Linearizing the viscous stress terms Tijj/p on the right-hand side of the momentum equation (la) and 
comparing magnitudes of the terms yields 


1 

P 



pU_ 


all other terms of being at least 0(h/L) smaller. 

In summary, the principal part of the viscous terms in the Navier-Stokes equations are 


iyP djjj + 2 pdijj'j 


( 7 - 1 ) 


pc 2 s, 


Pr \7 - 1 Cp 


- P,u 


( 7 - 1 )^ 
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pc 2 Pr \ 7 — 1 c p 


-P,i 


in the momentum, continuity, and entropy equations, respectively. 


(5a) 

(5b) 

(5c) 


3.3 The complete principal linearization 

Having found the principal linearization of the inviscid and viscous terms in Eqs. (4) and (5), the principal 
linearization of the full equations is found by comparing the relative magnitude of the terms in each equation. 

Starting with the entropy equation, a comparison of the magnitudes of the terms in Eqs. (4c) and (5c) 
yields the relative balance of advection and thermal diffusion 

c p M 2 ^ c p p U 2 
h pc 2 Pr h 2 

This can be restated as the requirement Re/,Pr = 0(1) for both terms to be principal. If the Prandtl number 
is 0(1), then the cell Reynolds number must be 0(1) for the viscous terms in the entropy equation to be 
principal. 

Examination of Eqs. (4a) and (5a) gives 


El >1— 

h p h 2 

as the dominant inertial and shear stress balance, or Re ft, = 0(1). 

Finally, the continuity equation terms from Eqs. (4b) and (5b) are comparable when 

2 U p U 2 
PC ~h ~ Pr ~h?' 

The relative magnitude of the viscous terms compared to the inviscid terms is M 2 /( Re^Pr), which for 
moderate Mach numbers corresponds to the condition that the cell Reynolds number is order unity. However, 
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for low Mach numbers the terms on the right-hand side of the continuity equation are subprincipal. Once 
again, this is reflection that in the incompressible limit, Eq. (lb) reduces to the statement that rqy = 0. 

To summarize, the principal linearization of the Navier-Stokes equations is 


1 . 

UjUij + -Pi 

p 

P C 2 Ujj + UjPj 


u jSj 





PC 


1) 


7-1 c p 

,2 


pc 2 Pr \7 — 1 
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P,u 


(6a) 

(6b) 
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3.4 Interpretation of the principal linearization 

As stated in the beginning of this section, the motivation for finding the principal linearization is that, 
because a relaxation method is effective at eliminating short wavelength error components, only those terms 
of the linearized operator that contribute strongly at short wavelengths are necessary. This property may be 
stated more formally. Let Lfl be the full linearization of the operator, and L be the principal linearization. 
The principal linearization satisfies the property 

HI-L^LflI \=o(^y (7) 

In other words, the inverse of the principal linearization is a good approximation to the inverse of the full 
linear operator at the scale of the grid spacing. 

An important point to emphasize here is that the key assumption used to derive the principal linearization 
is h <C L, which implies that the flow is well-resolved on the grid. It was shown above that formally the 
viscous terms are principal when Re^ = 0(1), which is quite restrictive. For example, a typical calculation 
of a laminar boundary layer flow over a flat plate will have on the order of 20-30 grid points across the 
boundary layer. The local Reynolds number in the middle of the boundary layer, based on the local flow 
speed and the normal grid spacing, is ~3— 10. When solving the system of equations using a full multigrid 
(FMG) process, the viscous region will be underresolved on the coarse grids. The principal linearization 
will not be strictly valid until the finest grid or two are reached. Also implicit in the condition h <C L is 
that the region of interest is well-removed from the boundaries and any discontinuities, such as shocks. The 
complications that arise near boundaries and discontinuies is discussed at length in Thomas et al. (2003). 

For typical aerodynamic flows where the viscous effects are confined to thin shear layers, the variation 
in length scales over the domain is quite large. Not all the terms will be principal in the different regions of 
the flow. In practice, the variation in physical scales is handled by variation in grid spacing throughout the 
computational domain. As a consequence, the principal linearization will differ throughout the computational 
domain. In the next section, various limiting cases of the principal linearization for Eq. (6) are considered. 


4 Factors of the Navier-Stokes equations 

In Section 3 the two components of an efficient multigrid method are given: a relaxation method for smooth- 
ing the short-wavelength component of the error, and a coarse-grid operator to reduce the long-wavelength 
error. For multigrid to be fully effective, it is necessary that the solution error is well-represented on the 
coarse grid. This is the case for elliptic partial differential equations. However, for hyperbolic equations not 
all smooth components of the error are visible to the coarse grid operator. Those components normal to 
the advection direction are only represented at the level of the truncation error. The coarse-grid correction 
for these contains only part of this smooth error component. Multigrid is not effective for these types of 
equations. 
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The physics of fluid flow is a combination of advection and diffusion phenomena, and as such the Navier- 
Stokes equations are a combination of elliptic and hyperbolic-parabolic operators. Energy and vorticity 
are carried along streamlines and transported transversely by diffusion. Conservation of mass represents a 
balance between pressure and the divergence of the velocity field, which is elliptic behavior subsonically and 
hyperbolic supersonically. To construct an efficient multigrid solver for the Navier-Stokes equations, it is 
necessary to split the system into these component parts such that each partition may be treated in the 
most appropriate way. 

The principal linearization in Eq. (6) can be used to identify and isolate the various subsystems of the 
Navier-Stokes equations. Introducing the notation Q = Uidi and Q u = Q — vA for convenience, the system 
may be written in the operator form L eq = r, or 


{Qv- {v + v')di\7- 

l di 
p 1 

0 ^ 

( euA 


pc 2 V- 

Q-^A 


e P 1 = r. 

(8) 

V o 

v c » {l - 1] A 

pc^rr 

Q~wA) 

w 



Because L is a matrix of operators, the formal eigenvalues of L are the “eigen-operators” of the system, in 
a generalization of characteristic analysis of hyperbolic equations. These eigenvalues are pseudo-differential 
operators. Unfortunately, not all these are true differential operators. However, taking the formal determi- 
nant of L gives an expression whose factors are identifiable differential operators. Some of these factors are 
the eigenvalues of L, while the remaining factors are products of those eigenvalues which do not correspond 
to any real differential operator. 

The determinant of Eq. (8) is 

det L = Qt 1 { (Q - ^ A) \q (Q - 7^ A) - c 2 a] +^-2u-^q(q~ 7 |- a) a} , (9) 

where d is the number of space dimensions. The repeated advection-diffusion factors correspond to the trans- 
port of vorticity in two and three dimensions, and to the conservation of helicity in three dimensions (Ta’asan, 
1993). The factor in the curly braces corresponds to energy and mass conservation. This factor is called 
the “core operator” by Brandt. Thermal diffusion couples temperature variations to pressure in Eq. (8). 
However, there are some special cases in which the core operator can be factorized, which will be examined 
below. Furthermore, some of these special cases will occur in regions of many flows of engineering interest. 
For this reason, these limits can be used when relaxing the full Navier-Stokes equations in these regions. 

4.1 Inviscid equations 

As shown in Section 3, the viscous terms are principal only where the cell Reynolds number is 0(1). For 
most aerodynamic flows, viscous effects are significant only in thin shear layers, and the bulk of the flow is 
inviscid. The inviscid limit of Eq. (8) is 

( Q o\ 

L„_> 0 = pc 2 V- Q 0 I ■ (10) 

V 0 0 QJ 

The advection-diffusion operator reduces to pure advection, and factors out of the core operator. The 
determinant of Eq. (10) is 

det L„^ 0 = Q d (Q 2 — c 2 A) . (11) 

Two of the advection factors correspond to vorticity and helicity conservation, as before. The additional 
advection factor is the transport of entropy along streamlines, which can be seen by inspection of Eq. (10). 
The remaining second-order operator is the full potential operator, which corresponds to conservation of 
mass. The full potential operator connects the pressure field to the irrotational part of the velocity field. 
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Note that in the one-dimensional case, the full potential operator factors into two first-order operators, 
corresponding to the left- and right-running acoustic characteristics. 

Although the discussion concerns the principal linearization of the equations, there is a correspondence of 
these terms in the full nonlinear case as well. When the flow is uniform at infinity, the vorticity remains zero 
and the entropy is constant. In such a case, all the advection terms are satisfied identically. The solution 
of the Euler equations is fully determined by the solution of the continuity equation alone. By writing the 
velocity as the gradient of a scalar potential function, the continuity equation reduces to the full potential 
operator acting on the potential function. Thus it is seen how the factors of the determinant of L in the 
inviscid limit correspond directly to a well-known property of inviscid, compressible flow. 

4.2 Low speed flow 

Another important limit is at low speed. Specifically, when M 2 <C 1, then the pressure advection term in 
the continuity equation is subprincipal. With fluctuations in entropy scaling as given in Section 3.1, the 
entropy diffusion term in the continuity equation is also subprincipal. In the entropy equation, the diffusion 
of pressure and entropy are of comparable magnitude, and remain principal. The resulting reduced form of 
the principal linearization is 


(Qv 




{v + 
pc 2 V- 
0 



—v 


pc 2 Pr 


A 



( 12 ) 


This corresponds to a nearly constant-density fluid, and occurs when the boundaries of the domain are 
adiabatic. 

If there is significant heating at the boundaries, then we can have that s, s c p M 2 . In this case, the 
temperature (or entropy) fluctuations are decoupled from the pressure fluctuations. If M 2 <C 1 the pressure 
advection-diffusion term in the continuity equation and the pressure diffusion term in the entropy equation 
are subprincipal, leading to the principal linearization 


(Qv-{v + v')di\7- ±<9i 


-*M— ► 0 — 
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pc 2 V- 

0 


A 

Cp Pr 

Q-& A 


For both Eqs. (12) and (13) the determinant of the operator is 


det Lm^o = -c 2 Q^ 1 (Q - — A) A. 


(13) 


(14) 


The decoupling of the temperature and pressure variations is reflected by the core operator factoring into 
an advection diffusion term and the Laplacian. The advection-diffusion factor acts on the entropy. The 
Laplacian correponds to conservation of mass. 

These two low speed limiting cases are distinct — one corresponds to incompressible flow, the other to 
low Mach number but compressible flow. An important point to note here is that the determinant of the 
operator is the same for both limits discussed above. By the criterion presented in Thomas et al. (2003), 
one would neglect altogether the off-diagonal dissipative terms in the continuity and entropy equations in 
the incompressible limit, resulting in the principal linearization 

/Q v -(v + i/)diV- jd t 0 \ 

Lm^o = pc 2 V- 0 0 

V 0 0 Q-&AJ 


This form of the operator does not satisfy the property given by Eq. (7). 



4.3 A special case of the Prandtl number 


Examination of Eq. (9) shows the term in the curly braces can be factored when the Prandtl number is 
1/(2 + p! jp). Substituting into the equations yields 


Lp r =l/(2 +h'/ij) 


(Q * 

V 


and this has the determinant 


(y + i /)diV- fa 

pc 2 X7 ■ Q~{l- l){2v + v')A 

0 -{2v + u') c -^^-A 


° \ 

~{2v + u')^ A 

Q-{2v + v')A J 


det L Pr=1/(2+M7 ^ = Qi 1 [Q - (2v + v') A] {Q [Q - 7 (2u + v') A] - c 2 A} . 


(15) 


(16) 


The new advection factor that appears reflects the transport of entropy along a streamline. The remaining 
term corresponds to continuity, but the diffusion terms reflect that the pressure is coupled to the entropy 
through thermal diffusion. 

In general, Stokes’ hypothesis is applicable which gives that p! = — 2/z/3, which means that this factor- 
ization applies when Pr = 3/4. Note that for air the Prandtl number ranges from 0.75 to 0.70 over the 
temperature range —100 to 100 Celsius, and is approximately 0.72 for most atmospheric flight conditions. 
Thus this factorization is very nearly exact for aerodynamic flows. 

For convenience, Eq. (15) and its determinant Eq. (16) is written below for the case when the Stokes 
hypothesis applies and Pr = 3/4: 


Lp r=3 /4 


(Q» - \vdiV- jd, 0 \ 

pc 2 V- Q — | (7 — 1)pA ~\^vA , 

V 0 A Q-¥N 


det L Pr=3/4 = Q d v 1 ( Q - jyA 


Q ( Q-lifA ) - c 2 A 


5 Modes of the principal linearization 


To construct an efficient relaxation scheme for the Navier-Stokes equation that makes effective use of the 
factors, it is necessary to determine the derived variables upon which the factors are acting. The most 
natural decomposition of the perturbation variables into modes corresponding to the factors of the operator 
is 




-p(Q - £ A ) 


v“(7- 1) 


c 2 Pr 


Ay 



(17) 


where U si is a solenoidal velocity field, and er and (j) are scalar potentials. This is related to the classical 
decomposition of a velocity field into solenoidal and irrotational parts. Applying the operator L to the above 
decomposition yields 


L 





7p? A ) T pc 2 A <j> 

0 / 



a 


+ (^-2,-*/)a 


( di 

V 


c D Pr 


0 

0 



(18) 


9 



Note that the last term vanishes for the special case of Subsection 4.3, i.e. , Pr = 1/(2 + p'/p). In this 
particular case, three remaining terms correspond exactly to the factors of Eq. (16). Furthermore, in the 
inviscid limit it is observed that cj) is the velocity potential, and a becomes the entropy. Most significantly, 
the continuity equation depends on alone, and the entropy equation depends on a alone. The coupling 
between these modes in the general case appears in the momentum equation, and arises from the effect of 
the dilatational viscous stresses. 

The quantity U si is an incompressible vorticity error transported along the mean flow. The quantity cr 
is an isobaric temperature error transported by the mean flow. The appearance of the velocity compo- 
nent v / (cpPr)di<r arises from the dilatation of the velocity field induced by heating of the fluid. The quan- 
tity </> is essentially a velocity potential. A change in cj) corresponds to a change in velocity at a constant 
total enthalpy, or total temperature. The continuity equation is the full potential equation modified by the 
addition of a term corresponding to the absorption of acoustic energy by thermal diffusion. In the case 
of Pr = 3/4, the velocity dilatation due to thermal diffusion is balanced exactly by the viscous action of the 
bulk viscosity. For any Pr ~ 1, the coupling is weak. 

In the inviscid limit, these three modes arc fully decoupled: 


0 [ U si\ 

= Q I 0 I 



4> + Q 



(19) 


It is interesting to note that these modes appear in Kovasznay’s study of supersonic turbulent flow (Clru & 
Kovasznay, 1958; Kovasznay, 1953). He developed the perturbation equations for small- amplitude, small- 
scale fluctuations being transported with a uniform mean velocity. In the present case, the short-wavelength 
solution errors are analogous to the small-scale turbulent fluctuations, and these errors are transported along 
the mean-flow approximate solution. Because the equations are linearized about an unconverged state, the 
residuals act as apparent mass sources, body forces, and heat sources on the right-hand side. Kovasznay’s 
work focuses on disturbances whose Reynolds number is much greater than one, whereas here the viscous 
terms are principal when Re/j = 0(1). For higher Reynolds numbers, the artificial viscosity of the scheme 
dominates the physical viscosity. 

It is useful to compare the present decomposition to that presented in Thomas et al. (2003). They use 
internal energy rather than entropy as their thermodynamic variable. Converting their formulation into the 
present notation, and rescaling some terms, yields 

V ) = ( p{v + v')V- ] U si + ( -p(Q - {2u + is') A) | <j> + ( 0 ) a. (20) 

s) V-T^ + ^V-/ \UQ-( to + v’)A)J \lj 


Here, the velocity U si is not assumed to be solenoidal. The quantity a is an isobaric temperature mode, but 
it does not induce any velocity dilatation. Finally, the quantity is not the usual velocity potential, but 
corresponds to a velocity and pressure change at a constant static temperature, rather than total temperature. 
The application of L to this set of variables yields 


L 



Q v \ ( o \ 

UV- ) U si + - pQ ( Q - (2v + u')A) + pc 2 A 

-^QV-J V ^Q{Q-{2v + v’)A) ) 


( ° \ 

( 21 ) 

\Q-wA) 


where V = pc 2 + (p + p')Q. Note that the energy and continuity equations are coupled for all values of Pr. 
Also, if the velocity U si is solenoidal, the tranport equation for this mode becomes fully decoupled from the 
other modes. In the inviscid limit, this system becomes 


fin A ( Q \ ( 0 \ (°\ 

L \ p = pc 2 V- U si + ~P (Q 2 - c 2 A) \ (j> + 0 cr. 

\s) \ 0 / V iQ 2 ) \QJ 


( 22 ) 
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6 Conclusions 


In this paper, the compressible Navier-Stokes equations have been examined with the objective of identifying 
the factors of the operator. The principle linearization of the systems has been developed by examining those 
terms of the linearization that are most significant for the short wavelength error components. The factors 
of the operator are seen to be identical to the modes identified by Kovasznay (1953) in his analysis of 
compressible turbulent flow. 

It must be emphasized that the decomposition of the Navier-Stokes equations into these factors is only 
the first step in the construction of a multigrid algorithm. A factorizable discretization of the system — one 
that preserves these factors at the discrete level — must be constructed. Such a discretization must also 
be /i-elliptic for an efficient relaxation scheme and multigrid solver to be constructed. Creating a scheme 
that preserves all these properties is an extremely difficult task (see Sidilkover, 1999). Indeed, there is no a 
priori guarantee that such a scheme exists. Nevertheless, the identification of the factors provides a useful 
framework for construction of fast converging schemes. 
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